N 6‐methyladenine profiling of low‐input multiplex clinical samples on transcriptome reveals RNA modifications implicated in type 2 diabetes and acute myocardial infarction

Dear Editor, RNA N6-methyladenines (m6A) not only play an essential role in normal biological processes but also participate in the pathogenesis of diseases including type 2 diabetes (T2D) and comorbidities.1,2 However, conventional technologies (e.g. m6A-seq)3 require large amounts of input materials often prohibitive for clinically-feasible biospecimens (e.g. a few millilitres of blood or tissue from fine needle aspiration), we, therefore, developed a Strategy of Low-Input Multi-barcode m6A-seq (SLIM-m6A-seq) to facilitate transcriptome-wide m6A profiling with limited RNA.Using this new technology, we implicated novelm6A modifications and pathways in the development of acute myocardial infarction (AMI) in patients with T2D. Technically, SLIM-m6A-seq exploited the barcoding strategy to encode multiple samples with unique barcodes before pooling, library construction, and sequencing (Figure 1A and Table S1).4,5 The resulting profiles can be split by barcodes (Table S2) and traced back to the original samples. We first tested the feasibility of SLIM-m6A-seq in HeLa cells and calibrated them with spike-in oligonucleotides (Note S1 and Figure S1). Then we compared SLIM-m6A-seq of a pooled sample containing barcoded mRNA varying from 200 to 10 ng with m6A-seq of 800 ng mRNA (Note S1). The SLIM-m6A-seq results from input levels of 200 and 100 ng showed similar distribution (e.g. chromosome 1 in Figure S2A) and complexity with 2.56% relative standard deviation (RSD), compared to the data fromm6A-seq in 800 ng input, albeit substantially different starting amounts (Figure S3A). We further validated the capability of SLIM-m6A-seq for low-input samples in two 25-sample sets containing 10 ng (Init10) or 5 ng (Init5) fragmented mRNA from HEK293T


N 6 -methyladenine profiling of low-input multiplex clinical samples on transcriptome reveals RNA modifications implicated in type 2 diabetes and acute myocardial infarction
Dear Editor, RNA N 6 -methyladenines (m 6 A) not only play an essential role in normal biological processes but also participate in the pathogenesis of diseases including type 2 diabetes (T2D) and comorbidities. 1,2 However, conventional technologies (e.g. m 6 A-seq) 3 require large amounts of input materials often prohibitive for clinically-feasible biospecimens (e.g. a few millilitres of blood or tissue from fine needle aspiration), we, therefore, developed a Strategy of Low-Input Multi-barcode m 6 A-seq (SLIM-m 6 A-seq) to facilitate transcriptome-wide m 6 A profiling with limited RNA. Using this new technology, we implicated novel m 6 A modifications and pathways in the development of acute myocardial infarction (AMI) in patients with T2D.
Technically, SLIM-m 6 A-seq exploited the barcoding strategy to encode multiple samples with unique barcodes before pooling, library construction, and sequencing ( Figure 1A and Table S1). 4,5 The resulting profiles can be split by barcodes (Table S2) and traced back to the original samples. We first tested the feasibility of SLIM-m 6 A-seq in HeLa cells and calibrated them with spike-in oligonucleotides (Note S1 and Figure S1). Then we compared SLIM-m 6 A-seq of a pooled sample containing barcoded mRNA varying from 200 to 10 ng with m 6 A-seq of 800 ng mRNA (Note S1). The SLIM-m 6 A-seq results from input levels of 200 and 100 ng showed similar distribution (e.g. chromosome 1 in Figure S2A) and complexity with 2.56% relative standard deviation (RSD), compared to the data from m 6 A-seq in 800 ng input, albeit substantially different starting amounts ( Figure S3A).
We further validated the capability of SLIM-m 6 A-seq for low-input samples in two 25-sample sets containing 10 ng (Init10) or 5 ng (Init5) fragmented mRNA from HEK293T cells. Overall, greater than 90% mapping rates on transcriptome were achieved for these samples (Table S3). A similar number of m 6 A peaks was detected, with an RSD of 34.26% in Init10 and an RSD of 37.43% in Init5 ( Figure S4A), as well as a consistent number of m 6 A-containing genes across different barcodes ( Figure S4B). The m 6 A profiles across the two sets were highly correlated ( Figure 1B and Figure  S4C,D), thus confirming technical reproducibility in lowinput samples. The barcoded samples in Init5 or Init10 sets exhibited high reproducibility of complexity with RSD ≦ 20% in both sets with or without m 6 A immunoprecipitation (IP), although the complexity would be reduced with smaller input as expected ( Figure S4E). The enrichment folds of m 6 A in individually barcoded samples were comparable between both sets -6.22-folds in Init5 and 7.51folds in Init10 ( Figure 1C and Figure S5A). The RRACH motifs ( Figure 1D and Figure S5B) and consistent m 6 A distributions ( Figure 1E) were discovered across all samples, as shown in an example of chromosome 1 ( Figure 1F).
We then profiled a cohort of 20 blood samples ( Figure  S6 and Table S4) from patients with T2D (n = 10) and T2D with AMI (T2D+AMI) (n = 10) matched for age, sex and blood pressure to demonstrate the clinical feasibility of SLIM-m 6 A-seq. A similar complexity was observed in all samples with the ratio of distinct/total reads ranging from 0.84 to 0.90 and an RSD of 1.74% ( Figure 2A). The principal component analysis suggested transcriptome-wide differences between these two groups ( Figure 2B). In total, 3223 dysregulated genes (|log 2 (fold change)| ≥ 1, adjusted p < 0.05) were detected between these two groups ( Figure 2C), with 1845 up-regulated and 1378 down-regulated in T2D+AMI ( Figure 2D). A panel of 14 differential genes was further established through the leave-one-out cross-validation and the elastic net regularization ( Figure 2E and Table S5). Note that, 7234 m 6 A peaks were found with typical distributions  ( Figure 2F) and RRACH motifs ( Figure 2G). SLIM-m 6 Aseq revealed an elevated overall m 6 A modification level in T2D+AMI patients ( Figure 2H), which was confirmed by high-performance liquid chromatography-tandem mass spectrometry ( Figure 2I and Table S6) and also reflected by the fact that 806 genes with increased and 260 genes with reduced m 6 A modification levels in T2D+AMI patients ( Figure 2J).
Notably, the expression level of the m 6 A demethylase FTO decreased evidently in T2D+AMI ( Figure 3A), supported by real-time quantitative polymerase chain reaction (RT-qPCR) ( Figure 3B and Tables S6 and S7). The conjoint analysis uncovered 392 genes exhibiting different combinations of expression and m 6 A modification patterns (|log 2 (fold change)| ≥ 1, adjusted p < 0.05 for expression change, and ≥ 1.68-fold enrichment for m 6 A). Of these intersecting genes, 311 genes exhibited concordant direction of changes, and 81 showed contrary alterations ( Figure 3C). Gene ontology (GO) analysis of the 392 overlapping genes revealed that the biological processes related to T2D+AMI were involved in protein kinase activities, inflammatory response, apoptosis, and angiogenesis ( Figure 3D), 6 which have been implicated in vascular biology, inflammation, and myocardial injury (Table S5). 7 Of the 392 intersecting differential genes, several inflammation-or atherosclerosis-related genes and angiogenesis activators showed correlation (Pearson |R| ≥ 0.5) with one or more clinical cardiac and inflammatory biomarkers ( Figure 3E). Of note, TNFSF14, a credible upregulated gene explored in our differential gene analysis and relevant to T2D and angiocardiopathy, 8,9 showed reduced m 6 A in T2D+AMI. This observation indicated that the effect of m 6 A in the progression of AMI during T2D could be through regulating the expression of disease-related genes. The GO analysis was further supported by the observed changes in clinical indexes relevant to myocardial injury and inflammation between T2D+AMI and T2D ( Figure 3F). Together with mediators in apoptosis and angiogenesis, these differential m 6 A-containing genes were relevant to inflammation and cardiomyocyte death, 10 suggesting the significant role of m 6 A in the occurrence of T2D+AMI through these mechanisms.
In summary, we developed SLIM-m 6 A-seq to implement transcriptome-wide m 6 A profiling with limited RNA. Comprehensive analysis of SLIM-m 6 A-seq data in clinical samples from patients with T2D and T2D+AMI provided novel insights into the epitranscriptome of T2D+AMI regarding the contribution of m 6 A. Although with certain limitations (e.g. potential confounding due to immune cell composition), the overall m 6 A levels in T2D+AMI appeared to be elevated compared to T2D, accompanied by down-regulation of FTO. Functional annotation analysis revealed relevant biological processes involved in the development of T2D+AMI. The current study demonstrated the technical robustness of SLIM-m 6 A-seq and importantly the feasibility of this new method in obtaining transcriptome-wide m 6 A from clinically convenient samples, thus opening up opportunities for exploiting m 6 A as an effective clinical approach for disease diagnosis and prognosis.

A C K N O W L E D G E M E N T S
We would like to thank experts in the department of cardiology (Professor Yanggan Wang, Hairong Wang and Jianlei Cao) and the department of endocrinology (Professor Haohua Deng) of Zhongnan Hospital of Wuhan University for disease diagnosis and critical discussion. The numerical calculations in this article were performed using the supercomputing system in the Supercomputing Center of Wuhan University.